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Abstract. The zero-temperature single-particle Green's function of correlated 
fermion models with moderately large Hilbert-space dimensions can be calculated by 
means of Krylov-space techniques. The conventional Lanczos approach consists of 
finding the ground state in a first step, followed by an approximation for the resolvent 
of the Hamiltonian in a second step. We analyze the character of this approximation 
and discuss a numerically exact variant of the Lanczos method which is formulated in 
the time domain. This method is extended to get the nonequihbrium single-particle 
Green's function defined on the Keldysh-Matsubara contour in the complex time 
plane which describes the system's non-perturbative response to a sudden parameter 
switch in the Hamiltonian. The proposed method will be important as an exact- 
diagonalization solver in the context of self-consistent or variational cluster-embedding 
schemes. For the recently developed nonequihbrium cluster-perturbation theory, we 
discuss the efficient implementation and demonstrate the feasibility of the Krylov-based 
solver. The dissipation of a strong local magnetic excitation into a non-interacting bath 
is considered as an example for applications. 
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1. Introduction 

Due to Wick's theorem, the single-particle Green's function G is the central quantity 
of interest in theoretical approaches to strongly correlated electron systems that are 
based on the concepts of weak-coupling perturbation theory [1, 2]. This holds for 
systems in thermal equilibrium as well as for systems subjected to strong time-dependent 
perturbations that give rise to a highly excited quantum state far from equilibrium 
[3, 4, 5]. Standard examples comprise the plain or different renormalized perturbation 
theories for lattice-fermion models like the Hubbard model [6, 7]. 
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Nonperturbative approximations can be constructed within dynamical variational 
approaches [8]. Here the single-particle Green's function or the self-energy is determined 
from a general variational principle. This includes dynamical mean-field theory (DMFT) 
[9, 10], different cluster extensions of the DMFT [11], the cluster-perturbation theory 
(CPT) [12, 13], the variational cluster approach (VGA) [14, 15] or the dual-fermion (DP) 
method [16], for example. These methods all involve a self-consistent or variational 
mapping onto an effective impurity or an effective cluster model for which G must 
be computed repeatedly. Among the standard "solvers" to treat those problems, such 
as the quantum Monte-Garlo method [17], for example, exact diagonalization or the 
Lanczos technique represents an important alternative. Due to the exponential growth 
of the Hilbert space with increasing number of degrees of freedom, the Lanczos solver is 
basically restricted to single-band models and comparatively small clusters or single-site 
approximations with a few orbitals per site only. 

As suggested by Gaffarel and Krauth [18], a two-step Lanczos procedure [19, 
20, 21, 22] can be used as an efficient method to get the zero-temperature Green's 
function G of the single-impurity Anderson model: After finding the approximate 
ground state of the model in a first Lanczos step, the frequency-dependent single- 
particle Green's function is obtained in the second step by approximating the resolvent 
{(jj—H)~^ I—)- '^jn{oj—Em)~^\'m){m\ with the help of the eigenenergies Em and eigenstates 
\m) of the Hamiltonian H in a, small Krylov subspace of the full Hilbert space. This 
Lanczos solver has turned out to be very efficient and reliable and is frequently used 
within DMFT, cluster DMFT and VGA, sec Rcfs. [9, 18, 23, 24] for examples. 

Its extension to the general nonequilibrium situation is, however, an open issue and 
represents the main motivation of the present paper. In detail our motivations are the 
following: 

First, we note that there is a growing need for theoretical nonequilibrium 
approaches to describe and understand recent experimental studies. This includes 
spin-relaxation and switching processes in nanostructured systems with itinerant and 
correlated electrons which are experimentally accessible by means of scanning-tunneling 
microscope techniques, for example [25, 26]. Another field of interest is given by 
fast-demagnetization processes probed by femtosecond optical excitations [27] or the 
nonequilibrium electronic structure of strongly correlated transition-metal oxides which 
may be monitored by femtosecond pump-probe spectroscopies [28, 29]. Furthermore, 
there is an urgent need for theory to understand the nonequilibrium dynamics of highly 
excited fermionic states realized in correlated systems of ultracold atoms in optical 
lattices [30]. 

Second, there are straightforward extensions of DMFT [31, 32] and its cluster 
variants, of the GPT [33] and the DF method [34] to correlated lattice models far from 
thermal equilibrium. Roughly speaking, these extensions are obtained when the theory 
is re-formulated on the Keldysh-Matsubara contour in the complex time plane. Different 
solvers for the resulting effective nonequilibrium impurity or cluster problem have been 
employed: An analytical approach based on the solution of a closed set of equations is 
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available in the case of the Falicov- Kimball model only [31]. For Hubbard-type systems, 
diagrammatic weak-coupling [32] and strong-coupling [35] perturbative techniques can 
be used, or variants of the continous-time quantum Monte-Carlo technique [17]. As 
exact-diagonalization solvers within nonequilibrium single-site or cluster methods, only 
full diagonalization procedures have been employed so far, namely for the nonequilibrium 
dual-fermion approach [34] and the nonequilibrium CPT [33]. Opposed to a Krylov- 
space construction, the Green's function G is calculated here from its spectral or 
Lehmann representation using a basis of the full Hilbert space of the effective impurity or 
cluster model. This limits the conveniently accessible system size to Lc = 6 sites/or bitals 
only which must be regarded as crucial as the convergence with Lc is known to be 
exponentially fast [18]. 

Third, in case of a sudden strong quench of a model parameter, i.e. far from thermal 
equilibrium, it is advisable to focus on real-time single-particle correlation functions 
of the form {ca{t)cl^/{t')). Here a, a' refer to one-electron orbitals and (• • •) is the 
expectation value with an initial state different from the ground state or an eigenstate 
of the time-independent Hamiltonian H. In this case, the question for an approximation 
of the operator exponential exp{—iHt) rather than the resolvent must be addressed. A 
standard and rehable Krylov-space method is available to propagate a given state |^) 
via exp{—iHt)\'^) [36, 37, 38, 39, 40]. The method is also used in the context of the 
density- matrix renormalization [40] and continuous-time quantum Monte-Carlo [41]. In 
the present paper, we focus on an efficient application of this Krylov approach to get 
the Green's function on the Keldysh-Matsubara contour with a maximum real time 
^max and a time discretization step At typical for solver applications. To demonstrate 
its feasibility, we employ the approach within the context of nonequilibrium GPT. 

Finally, for the equilibrium single-particle Green's function G, the Krylov 
construction to approximate the exponential exp{—iHt) represents an alternative 
approach to the conventional Lanczos technique which approximates the resolvent 
(u — H)~^. Their mutual relation shall be worked out here. 

The paper is organized as follows: The next section gives a brief overview of the 
Krylov construction and the standard Lanczos approach for the equilibrium Green's 
function. The approximate nature of the approach is made clear and the type of 
the approximation is characterized. Section 3 discusses a variant of the technique 
formulated in the time domain by which the numerically exact Green's function is 
accessible. In section 4 we seek an algorithm to get the nonequilibrium Green's function 
on the Keldysh-Matsubara contour in the complex time plane and propose a four-step 
Krylov-space based technique. To demonstrate its feasibility in the context of a cluster- 
embedding scheme, we consider the dissipation of a local magnetic excitation in section 
5 by means of nonequilibrium CPT. Section 6 summarizes the main results. 
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2. Krylov space basis and equilibrium Green's function 

We start by giving a brief overview of the Lanczos approach to the single-particle Green's 
function. Details can be found in Ref . [22] , for example. 

For a given initial state |io), the n-th Krylov space is defined as 

)Cn{\to)) = spe.n{\to),H\to),...,H"-%)}, (1) 

where H is the Hamiltonian of the system. /C„(|io)) is an n-dimensional subspace of the 
full Hilbert space with dimension d. Usually, we consider n <^ d. A basis of ICn{\io)) 
can be constructed by means of the numerically efficient recursion scheme 

\ik+i) ^ H\ik) - ak\ik) - bl\ik-i) (/c = 0, n - 1) (2) 

with initial values Bq = and = and with the coefficients = {ik\H\ik) / (ikl'ik) 
and bl = {ik\ik) / {ik-i\ik-i) ■ Subsequent normalization yields the orthonormal Lanczos 
basis {|io), Kn-i)} of /Cn(Ko))- In this basis, the Hamiltonian is represented 
by a tridiagonal matrix T with diagonal elements given by Q-n-i and off- 

diagonal elements by 6i,...,6„_i: Let H be the d x d matrix representation of the 
Hamiltonian in an arbitrary basis {\j)}, e.g. in the occupation-number basis where 
\j) — |ni, 712, tIq, ...) and Ua are the occupations of single-particle orbitals |q;). We 
have Hjji — {j\H\j'). Let V = (io, in-i) be the d x n matrix constructed from 
columns ik representing \ik) in the given basis, i.e. Vjk — {j\ik)- Then 

T = V^HV . (3) 

Diagonalization of T, 

D = Q^TQ (4) 

with a unitary n x n matrix Q, yields a diagonal matrix D containing approximate 
eigenenergies Em of H. The corresponding approximate eigenvectors, i.e. H\m) ~ 
Ejn\m) for m = 1, n, are 

\m)^Y.^,^JJ) (5) 

i 

where we have defined the dx n matrix U = VQ. 

The convergence of the extremal eigenenergies with increasing n is very fast. To 
get the ground-state energy and the ground state itself, numerically almost exact results 
can be obtained with of the order of n = 100 Lanczos iterations [21]. The initial state 
\io) is arbitrary but must have a finite overlap with the ground state. 

Consider now the single-particle Green's function. For frequency > 0, the zero- 
temperature retarded Green's function is given by 

GMu^) = Gi>J,{u) = {0\c ^ff , F ^I'IO) (6) 

U + 17] — si + -t/Q 

where |0) is the ground state which is assumed to be nondegenerate, Eq is the ground- 
state energy, Cq, annihilates a fermion in the one-particle orbital la), and is a small 
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positive number to shift the poles of the Green's function below the real axis in the 
complex frequency plane. Forcj < 0, Gaa'i^) = G[^J,{uj) = {0\cl^,{u+iri+H—EQ)~^Ca\0)- 
The Lanczos procedure [21, 22] to get the Green's function consists in the following 
approximation for the resolvent (for cu > 0): 

77 7^ 7^\'m') ('ITT'l (7) 

m 

where the n <^ d approximate energy cigcnstatcs \m) arc obtained with Eq. (5) from a 
second Lanczos run using |io) = c]!^'|0) as the initial state. Typically, n ~ 100 is used 
again. This yields a Lanczos Green's function with exactly n poles: 

G^A^) ^ Gi^^,{u) ^ Y.(^\ca\m) \ H4|0) . (8) 

^ uj + ir] — him + -C/O 

A continuous spectral function A^a'ii^) = — (l/7'')IiriGQ,Q,/(a;) is obtained by Lorentzian 
broadening with a finite 77 > O. 

To estimate the quality of the approximation (7), we consider the high-frequency 

expansion of the exact Green's function (6) for 77 — > O, 

00 ^ 

G^^u) = J2 ^,{^K{H - EoYcllO) , (9) 

r=0 

and compare with the high-frequency expansion of the Lanczos Green's function (8). 
Introducing P = ^ projector onto /C„(|io)) with \io) — c^/|0), we 

immediately get for the latter: 

00 ^ 

g2h = E - Eo)prci\o) . (10) 

r=0 

Here, we have assumed that the error in the determination of the ground state |0) can be 
neglected. This is usually an excellent approximation which will also be adopted in the 
rest of the paper. Comparing Eqs. (9) and (10) shows that the Lanczos approximation 
for the Green's function conserves the first n coefficients in the high-frequency expansion 
since {H - EoYcl\0) = {H - EoY\to) = {P{H - Eo)PY\to) e }Cn{\to)) if r < n - 1. 

The expansion coefficients determine the first n moments J^^dix! uj^Aaa'i'-^) of the 
spectral function. Therefore, we can conclude that the Lanczos technique at iteration 
depth n provides a spectral function with the correct first n moments — irrespective 
of the fact that the excited states |m) are obtained with a much lower accuracy than 
the ground state. However, significant deviations from the exact spectral function are 
expected at high excitation energies since the convergence with increasing n is known 
to be faster for low-lying as compared to highly excited states. 

3. Numerically exact computation of the Green's function 

Seeking for an improved approximation, let us consider the time-dependent Green's 
function 

1 f°° 

Gaa'it) - — / du e-'^^'G^^iiu) (11) 
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which is obtained from the frequency-dependent Green's function in Eq. (6) 
via Fourier transformation. A straightforward calculation yields Gaa'{t) — 
—iQ(t) (0|cae^'*^^^'^"''*cj^, |0)e^''*. At this point we can employ a Krylov-space technique 
[36, 37, 38, 39, 40] to compute the time evolution of the state 

|*«,(i))^e-^H,|0), (12) 

and therewith 

G^ait) = -ie{t) (0|c,|vl/„,(t))e^(^°+^^)* . (13) 

The idea is the following: If At is sufficiently small, the Taylor expansion of 
the operator exponential exp{—iHAt) can be truncated at some finite small order 
n <^ d within numerical accuracy. This implies that the state \^a'{t + At)) = 
exp(—iHAt)\^a'{t)) lies in the Krylov space }Cn{t) constructed from the initial state 
|io) = |^Q,'(t)) at time t (for simplicity we suppress the a' dependence of /C„(i) in the 
notation). Hence, 

\-^a'{t + At)) = eM-iPit)HP{t)At)\-^a'{t)) , (14) 

where P{t) is the projector onto /C„(t). Within this Krylov space, the time evolution 
operator can be represented as 

exp{-tP{t)HP{t)At) = J2 \j)[exp{-tV{t)Tit)V\t)At)y{f\ , (15) 

jj' 

where V{t) and T{t) are the representations of the Lanczos basis and the Hamiltonian 
obtained during the construction of /C„(i) via the Lanczos iteration. Using the 
orthonormality of the Lanczos basis, V^V — 1 VV^, and Eq. (4), we have 

eM-iPmPit)At) = J2 \j)P{t) exp{-iDit)At)U\t)]jy{j'\ , (16) 

jj' 

where D{t) = Q{tyT{t)Q{t) and U{t) = V{t)Q{t). 

For sufficiently short At, Eq. (16) thus provides a numerically exact way to 
propagate the state (12) by At using the Lanczos recursion algorithm. The time 
propagation can be repeated by restarting the algorithm with the state at t + At as the 
new initial state. Using several restarts, this allows us to compute the time-dependent 
Green's function from t = up to a time tmax at which, depending on the choice of rj, 
the exponential damping e"''* ensures convergence of the time integral in the Fourier 
back transformation. This yields the frequency-dependent Green's function G^aK'^) 
for a given rj. It is worth mentioning that this approach ("time-dependent Lanczos") 
is numerically exact, opposed to the "conventional" Lanczos procedure described in 
the preceding section. The essential difference is that the computation of the operator 
exponential exp{—iHt) can be decomposed into several steps with short At, and that 
exp{—iHAt) can be represented numerically exactly by means of a low-dimensional 
Krylov space. On the other hand, the resolvent l/{uj — H) must be computed in a single 
step and thus be represented in a single (larger) Krylov space. 
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Figure 1. Local spectral density A{uj) — {—l/Tr)lmGu{uj) of the particle-hole 
symmetric Hubbard model on an i = 10 site ring at C/ = 8. The nearest-neighbor 
hopping T — 1 sets the energy scale. Lorentzian broadening: rj ~ 0.01, maximum time 
for propagation of states: imax — 1000.0. Inset: numerically exact spectral function 
obtained with / = 1000 Lanczos restarts, i.e. At — 1.0. Left panel: f = 1, At — 1000.0 
(red fat line) compared with / = 1000, At ^ 1.0 (blue thin line). Right panel: f = 100, 
At = 10.0 (green fat line) compared with / — 1000, At — 1.0 (blue thin line). 



Figure 1 gives an example for the one-dimensional Hubbard model 

H = -TJ2Y1 + UJ2 ni^n,^ (17) 

(ii> o-=t,i » 

with nearest-neighbor hopping T = 1 and Hubbard interaction f/ = 8 at half-filling. 
We consider a system with L = 10 sites and periodic boundary conditions. The local 
retarded Green's function Gii{uj) is calculated via numerical fast Fourier transformation 
from the time- dependent Green's function Gii{t). The latter is obtained via Eq. (13) for 
times up to tmax = 1000.0 which is sufficient to ensure the convergence of the Fourier 
transform at a Lorentzian broadening of rj = 0.01. 

For the propagation of a state by means of the Krylov technique, Eq. (12), it is more 
advantageous to consider a large Krylov space dimension n and a large time propagation 
step At but a smaller number of restarts / opposed to a small n and short At but more 
restarts /. On the other hand, the Krylov space dimension should not be much larger 
than n = (9(100) since all states of the Lanczos basis (i.e. V) have to be stored. For 
the present calculation, we fix n = 100. To avoid a loss of orthogonality of the Lanczos 
basis states during the iterative procedure, a Gram-Schmidt reorthogonalization scheme 
is employed. 
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The spectrum obtained with / = 1000 restarts, corresponding to At — 1.0, is shown 
in the inset of figure 1. This represents the numcricaUy exact solution, any increase of 
/ or n does not change the results. Note that Eq. (16) can be used to compute the 
Green's function on a finer time grid with spacing At' < At without additional restarts. 
Here we have used At' = 0.01 independent of the different / considered. In the left 
panel of the figure the numerically exact result is compared with the spectrum obtained 
for / = 1 (due to particle-hole symmetry only frequencies a; > are displayed). Here, 
the Krylov space is constructed only once, i.e. At — tmax = 1000.0. As can be seen from 
the figure, there is a perfect agreement for lower frequencies while deviations are clearly 
visible for the number and the energy position of the peaks as well as for their spectral 
weights at higher frequencies. 

It is important to realize that the / = 1 spectrum just corresponds to the result 
of the conventional Lanczos method: If the Krylov space is constructed only once from 
the initial state \io) = c]^/|0), the matrices D and U in Eq. (16) are independent of t. 
Inserting Eq. (16) with At replaced by t into Eq. (14) and Eq. (13), yields 

Gaa'it) = -^e(^) 5^(0|c„|m)e-*^-*(m|cl,|0)e'(^°+*^)* , (18) 

m 

where m just runs over the approximate eigenstates, see Eq. (5), that are obtained from 
the conventional Lanczos technique. Fourier transformation of Eq. (18) gives G^^^,{ijo) 
as defined in Eq. (8). We arrive at the conclusion that time-dependent Lanczos carried 
out with a single Krylov space (/ = 1, see figure 1) is equivalent with the conventional 
Lanczos method. This has also been checked numerically. 

The left panel of figure 1 therefore shows the deviations of the conventional Lanczos 
method from the exact result. The perfect agreement at low frequencies is now easily 
explained by the fact that the ground state and the low-lying excited eigenstates of H 
are accurately predicted by the conventional Lanczos method. Discrepancies at higher 
frequencies of the order of U are attributed to the poor convergence of higher excited 
states. At even higher excitation energies, outside the support of the spectrum, the 
conventional Lanczos Green's function becomes reliable again since the first n moments 
and thus the corresponding coefficients in the high-frequency expansion are predicted 
correctly as discussed in section 2. Note that the time-dependent Lanczos approach 
does not make any reference to the excited eigenstates of H although for a large Krylov- 
space dimension, such as n = 100, the elements of D{t) and of U{t) may be close to 
the eigenenergies and to the coefficients of the eigenstates [Eq. (5)] and only weakly 
dependent on the initial state for the Lanczos restart at time t. On the other hand, this 
weak dependence on the initial state is important to get the numerically exact result. 

The time-dependent Lanczos is as memory efficient as the conventional one. Since 
all states of the Lanczos basis have to be stored (the matrix V), memory requirements 
are minimized for a small n. Very small Krylov-space dimensions (e.g. n < 10) may be 
used at the cost of an increased number of restarts / (i.e. short At). On the other hand, 
CPU time is minimized with a small / and large n. As compared to the conventional 
method, the computational cost is to a very good approximation higher by a factor / 
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(for the same n) since the Krylov space must be constructed / times. This must be 
kept in mind for apphcations hke DMFT. 

This raises the question for a possible compromise: Can a small number of restarts 
/ cure the errors of the conventional Lanczos approach? The right panel of figure 1 
shows the result of the spectral function from a calculation with only / = 100 restarts 
keeping the Krylov-space dimension unchanged (n = 100). It turns out that the time 
evolution of the state, Eq. (14), is no longer exact over the entire time interval At — 10.0 
for each restart. Compared to the conventional Lanczos method, the deviations from the 
exact spectral function are of the same order of magnitude. More important, however, 
the approximation is no longer causal and produces negative spectral weight as can be 
seen from the figure. The case / = 1 represents an exception. Here, Eq. (8) applies and 
the non-negativity of the local spectral function is obvious. 

4. Nonequilibrium Green's function 

The nonequilibrium single-particle Green's function depends on two time arguments and 
is given by 

Gaa'iz, Z') = -i{0\rCaiz)cliz')\0) (19) 

where |0) is an arbitrary state, usually not an eigenstate of H, which describes the 
system at a time t — 0, and where the annihilator and the creator are given in the 
Heisenberg picture with times z,z' on the Keldysh-Matsubara contour in the complex 
time plane [5]. T denotes the time ordering on the contour. We are seeking an algorithm 
to compute the Green's function by means of a Krylov-space technique that meets 
the requirements for a "solver" in the context of nonequilibrium DMFT [31, 42] or 
cluster-embedding approximations, such as the nonequilibrium CPT [33]. This means 
that impurity or cluster models at half-filling with more than L = 6 sites should be 
accessible. In any application as a solver, Dyson's equation, which is an integral equation 
on the contour, must be solved by time discretization and numerical matrix inversion. 
Therefore, the Green's function must be computed on a discrete time mesh on the 
contour that is sufficiently fine for applications of standard quadrature formulas. Finally, 
since inversions of matrices in the time variables are involved, the typical maximum time 
^max up to which the time propagation of observables is traced is comparatively small, 
e.g. tmax = 10 in units of the inverse hopping. 

Let us briefly discuss the possibility to compute the Green's function in frequency 
space where, like in the conventional Lanczos approach, the resolvent is approximated. 
Consider, for example, the lesser Green's function 

G^^M = z{0\cl,{t')c^{t)\0) = z(0[e^^*'cl,e'^(*-*')c„e-^*[0) (20) 

with real time arguments t,t' > and t — t' > 0. Using the identity J^^duje~^'^^{uj + 
ir] — H)~^ — —2mQ{t)e~^^^e~^^, this can be written as 

Gta'M = j j j dw,duj2du;s e--it'g-ia..(t-t')e--3t ^ 
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X (01 ^ c^ ^- c„ ^- 10) . (21) 

Now, using |0) as the initial state for the Lanczos iterations, one may construct the 
Krylov space /C„(|0)) and, using Eq. (5), the orthonormal basis {|?Ti)} of /Cn(|0)) 
consisting of approximate eigenstates of H. Therewith the resolvents (cus + irj — H)^^ 
and {oji + it] + H)~^ can be approximated like in Eq. (7). For the remaining resolvent 
{uj2 + ir) + H)~^ , another basis must be constructed for any m if, in the spirit of the 
conventional Lanczos approach, Ca\m) shall be used as the respective initial state. Even 
then, however, the high-frequency asymptotics cannot be recovered correctly, opposed 
to the equilibrium case. This final step, therefore, represents a crude approximation. 

A numerically exact access to |\I'(a;3)) = {u-^ + it] — H)^^\0) (and likewise to 
(0|(a;i + irj + H)~^) would be provided by the correction-vector method [43, 44] which 
is frequently employed in the context of dynamical density-matrix renormalization 
[45, 46, 47, 48]. For each frequency cua, the correction vector \^{u!s)) can be obtained as 
the solution of a sparse inhomogeneous system of linear equations [cus+ir] — H)\'i/ (cus)) = 
|0) with a dimension given by the Hilbert-space dimension. To evaluate Eq. (21), 
however, another correction vector, depending on two frequency arguments, must be 
computed as the solution of {uj2 + if] + H)\x{uj2, (^3)) = Ca|^(a;3)). This appears as less 
efficient than approaches working in the time domain directly. 

We therefore propose the following four-step procedure to compute the t, t'- 
dependent Green's function: 

(i) The system's initial state |0) must be given or is calculated by means of 
a standard Lanczos procedure as the ground state of an initial-state Hamiltonian 

-f^ini 7^ H. 

(ii) Constructing the Krylov space with |0) as the initial state for the Lanczos 
iteration and, depending on n and imax, using / additional restarts, the state 

mt)) = e-^^*|0) (22) 

is computed with numerical accuracy and stored on a discrete time mesh for all times 

up to tmax- 

(iii) For each orbital a of interest, the state 

|M/,(t))^e^^*cJ$(t)) (23) 

is computed with numerical accuracy and stored on the time mesh for all t up to tmax- 
This step is most time consuming since the Krylov time evolution must be performed 
for any t on the time mesh, i.e. for any initial state Ca\^{t)). The CPU time for the 
construction of a Krylov space scales linearly with the dimension n. While usually it 
is efficient to employ large Krylov spaces and longer propagation times At, one has to 
bear in mind that the Green's function (19) must be obtained on a fine time mesh which 
requires, see Eq. (16), the computation of 0{n) dot products of Hilbert-space vectors 
for each time on the mesh. Therefore, it is not advisable to construct too large Krylov 
spaces. With increasing n, the CPU time is eventually dominated by the evaluation of 
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Figure 2. Lesser and greater local spin-t Green's function Gii{0,t) and Gfi(t,0) as 
functions of t for a single- impurity particle-hole symmetric Anderson model with L = 8 
sites in a chain geometry with the first site (i = 1) as the correlated impurity and with 
U = 1. Energy and time units are set by the nearest- neighbor hopping T = 1. At 
t ~ the system is prepared in a state that is given by the ground state of the same 
model but in the presence of a local magnetic field of strength h — 0.2 applied to the 
impurity site. The field is switched off for t > 0. 



Eq. (16) rather than by the Krylov-space construction for the case that there are many 
time points within a single interval At. 

(iv) Finally, the lesser Green's function at arbitrary times t, t' is obtained as the 
scalar product: 

G<At.t')=i{^ Amo.it)). (24) 

Likewise G^^,{t,t') and Green's functions G\^,{t,T') and G^^^,{T,t') with mixed 
real/imaginary time arguments can be calculated while the Matsubara Green's function 
G^„,(r, r') = G^^/(r — r') is the equilibrium Green's function for the initial-state 
Hamiltonian ifjni and accessible by the time- dependent or by the conventional Lanczos 
method described in section 4. 

Figure 2 gives an example for the single-impurity particle-hole symmetric Anderson 
model on a one-dimensional chain with L = 8 sites and with the correlated impurity at 
site 1: 

L-l 

H = Simp ^I'^^l^ + ^ '^It^U - ^ 5Z 5Z {^I'rCi+la + H.C.) . (25) 
a a i=l a 

Here, T is the nearest-neighbor hopping, U the Hubbard interaction, eimp = —U/2, and 
nia = cjo-Cio-, 0" i- The system's state |0) at time t = is defined as the ground state 
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of an initial-state Hamiltonian Hi^^ = H + i^geid which includes a finite local magnetic 
field at the impurity site in addition. The field disturbs the system at t — and is 
switched off for t > 0: 

^^fieid(^) = -he{-t)(ni^ - nij . (26) 

The initial state |0) is calculated by means of the standard Lanczos technique with 
n = 200. For t = 0, we have ImG< (0,0) = (0|c5^cit|0) and ImG'fi(0,0) = -{0\ci-ic\^\0) 
for the spin-t lesser and greater Green's functions, respectively. Their difference is unity 
as can be seen in the figure and as required by the canonical anticommutator relations. 
The real parts must vanish. For t > 0, Gfi(0, t) and 0) become complex and show 

strong oscillations as it is typical for a finite-size system. 

We have compared the suggested Krylov-space method to compute the Green's 
function (19) on the complete Keldysh-Matsubara contour with a fuU-diagonalization 
approach. With the latter all eigenstates of the Hamiltonian are obtained by numerical 
diagonalization. Time dependencies and Gaa'iz, z') are easily obtained then. Typically, 
the fuU-diagonalization approach is faster up to L = 6 sites at half-fiUing and exploiting 
the symmetries due to conservation of the total particle number and the z-component 
of the total spin. For L = 8 and larger systems, the Krylov approach is superior. 
Eventually, both the full diagonalization and the Krylov method are limited by the 
need to store matrices of size d x d or d x n, respectively. At half-filling L = 12 sites are 
easily accessible with the Krylov method on a standard PC. CPU times are an order of 
magnitude longer as compared to the standard Lanczos technique for the equilibrium 
Green's function. 



5. Nonequilibrium cluster-perturbation theory 

The nonequilibrium cluster-perturbation theory (NE-CPT) [33] is a simple cluster- 
embedding approach and constructed as a straightforward generalization of the standard 
(equilibrium) CPT [12, 13]. The NE-CPT can in principle be applied to an arbitrary 
lattice model of correlated electrons with local interactions. The main idea is to partition 
the original lattice into smaller pieces (clusters) for which the Green's function can be 
computed by means of an exact-diagonalization approach. The Green's function G 
of the original model is then obtained from the cluster Green's function G', or more 
precisely the Green's function G' for the system of disconnected clusters, via the CPT 
equation: 

G = ^ . (27) 

G'-^ -V ^ ^ 

Here, V is the mter-cluster hopping. The CPT equation can be interpreted as 

a resummation of diagrams in a perturbative expansion of G around the limit of 

disconnected and non- interacting clusters [33]. Thereby certain vertex corrections are 

neglected which describe the effects of inter-cluster potential scattering on the electron 

self-energy. In the diagrammatic formulation, the step from the CPT to the NE-CPT is 
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t=t 
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• — o — o — o- - -o — o — 

-hT T T T 



infinite batli 



t>t 







U V 

•— o — o — o- - -o— o — 

i=iT 2 T T Lc T 



Figure 3. Pictorial representation of the final-state Hamiltonian H governing the 
time evolution (bottom) and the Hamiltonian Hi^y^ generating the initial state as its 
ground state (top). Note that the Hubbard interaction U is spatially separated from 
the hopping perturbation V which links a small cluster of Lc < 12 sites with an infinite 
uncorrelated bath in a semi-infinite chain geometry. The initial state is defined to be 
the ground state in the presence of a finite local magnetic field h — 0.2 applied to 
the impurity site. The field is suddenly switched off at <o — 0. We consider the 
particle-hole symmetric model at half-filling. 



particularly clear since standard perturbation theory for a nonequilibrium situation 
basically follows just along the lines of perturbation theory for systems in thermal 
equilibrium [5]. As concerns the CPT, the essential new point is that all quantities in 
CPT equation (27) have to be interpreted as given on the Keldysh-Matsubara contour in 
the complex time plane, see equation (19), and that the matrix inverse in (27) not only 
refers to the orbital indices but also to the time variables, i.e. G is actually given by the 
solution of an integral equation. While the conventional Lanczos method is frequently 
used as a solver for equilibrium CPT [12, 49, 50, 51, 11], only full diagonalization has 
been considered for NE-CPT so far. Half-filled Hubbard clusters of no more than Lc = 6 
sites can thereby be treated conveniently. 

Usually, the disregard of the mentioned vertex corrections represents a severe 
cluster mean-field-type approximation. The correction to the self-energy of lowest 
order in V , however, depends on the square and higher powers of the free off- diagonal 
Green's function that links sites with finite U and V interactions. In cases where these 
interactions are spatially separated, vertex corrections are expected to be small. We 
therefore consider a single- impurity Anderson model (SIAM) on a semi-infinite chain 
with the impurity on the first site but with the "inter-cluster" hopping V between sites 
Lc and L^ + l (see section 4 and figure 3). V connects a small SIAM with < 12 sites 
and an infinite uncorrelated bath. The Green's function for the disconnected system 
with V = therefore consists of two independent parts: the Green's function of the 
isolated cluster and the bath Green's function. 

The bath Green's function G^^'^ is readily obtained as 




(28) 



if z later than z' on the contour and 




(29) 
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Figure 4. Left: Local magnetic moments rrii = Ui^ — in tlie initial state at t = 
and as functions of time t > at the cluster sites i = 1, ...,Lc = 8 for the decoupled 
{V — 0) model displayed in figure 3. Right: Profile of the moment at i = 10. Inset: 
mi for i > 5 on a larger scale. Results are obtained from the nonequilibrium Green's 
function, riia-it) = —iGfi^{t,t), which has been calculated using the four-step Krylov- 
space method. Parameters: Lc — 8, T = V — U — I, h = 0.2 for t = 0, h — ior t > 0. 
The initial state has been obtained as the ground state of 7?ini using the conventional 
Lanczos technique with n = 200 and using Gram-Schmidt reorthogonalization. For the 
final-state dynamics At ~ tmax = 10 (no restart) and different Krylov-space dimensions 
(as indicated) have been used. 



if z' later than z. Here, the parameter /3 — > oo projects out the ground state of the bath 
Hamiltonian with hopping matrix Tb. 

The Green's function of the isolated cluster is calculated using the Krylov-space 
method discussed in the preceding section. Its dependence on the time variables is 
inhomogeneous and reflects the system's time evolution after an initial perturbation. 
The initial state is taken to be the ground state of the system but with a finite local 
magnetic field apphed to the impurity site with strength h = 0.2, see equation (26). 
This field polarizes the vicinity of the impurity site in the finite cluster. Figure 4 gives 
an example for a cluster size Lc = 8. At t = there is a strong local magnetic moment 
at the impurity site. With increasing distance from the impurity the moments alternate 
around zero and decrease in size. Note that the total polarization Yli=i = as the 
field is too weak to break up the singlet ground state of the cluster. 
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Figure 5. Time dependence of the polarization of the first four sites for the model 
displayed in figure 3. Results for isolated {V — 0) clusters of different size Lc (dashed 
lines) and for the full {V = T) system (solid lines) as obtained by NE-CPT for 
Lc = 4 - 10 at t/ = T = 1. 

The field is suddenly switched off at t = 0. The site-dependent moments for times 
t > are obtained from the time-diagonal elements of the cluster Green's function 
(opposed to the off-diagonal elements shown in figure 2). As expected physically, the 
strong impurity polarization dissipates into the rest of the system and decreases with 
time. However, the system is finite and small which causes a strong revival of the 
moment at a time ~ 10. For t ^ 5 one can see the polarization to be at a maximum at 
2 = 8, i.e. at the opposite edge of the chain. 

A maximum time of the order of tmax = 10 is dictated by the CPT due to the 
necessity to solve the CPT equation by time discretization and inversion of matrices 
in t,t' (see also below). This limitation of the CPT to the short-time physics is 
actually characteristic for any nonequilibrium cluster-embedding method although a 
somewhat larger tmax is possible using advanced quadrature formulas. For tmax = 10 
the calculations can be done by constructing the Krylov space only once without any 
restart. Namely, as can be seen by comparing the results for different Krylov-space 
dimensions n in figure 4, convergence is obtained for 50 in this example, i.e. for 
rather moderate values. 

With the contour-ordered Green's function of the isolated cluster and the bath 
Green's function at hand, i.e. with G', the Green's function of the full model G can 
be obtained from (27). According results are shown in Figure 5. For the calculation, 
we have considered the SIAM with a cluster of = 4 — 12 sites and a bath consisting 
of Lb = 1000 sites simulating a semi- infinite system (see figure 3). The full contour- 
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ordered Green's ftmction Gij„{t,t') for sites i,j in the cluster is obtained by solving the 
time-discrctizcd CPT equation (27). For the Kcldysh branch along the real time axis 
from t = up to t = tmax = 10, a time integration step of 6t = 0.02 has turned out to 
be sufficient for convergence. For the Matsubara branch along the imaginary axis we 
find converged results for /3 = 15 and 5t = 0.02. Note that not only the time evolution 
of the final state but also the initial equihbrium state is treated by means of the CPT. 
Therefore, to take potential-scattering vertices into account for the initial state, the 
Matsubara branch must be included in the calculation. For further technical details on 
the NE-CPT and the solution of the CPT equation exploiting symmetries, we refer to 
[33]. 

For the discussion of the results, consider the initial state ai t — first. The 
magnetic field h = 0.2 applied at the impurity site causes a sizable impurity polarization 
mi = riif — rill ~ 0-24 (not shown). For the half-filled system, strong antiferromagnetic 
spin correlations then lead to a polarization cloud close to the impurity with alternating 
local magnetic moments rrii which decrease in size with increasing i. The figure shows 
the net polarization at the first four sites of the semi-infinite system which amounts to 
Yli=i "^j ~ 0.137 in the initial state. Note that the convergence with increasing cluster 
size Lc is extremely fast for the CPT results (full blue fines at i = 0) as compared to 
the results for an isolated cluster (dashed orange lines ait — O). For the isolated cluster 
with Lc = 4 the total polarization even vanishes as the ground state is a singlet for the 
considered small h. 

The finite net polarization close to the impurity is expected to dissipate into the 
infinite, initially unpolarized bath. Indeed, the CPT shows that as a function of time the 
polarization relaxes quickly, overshoots a bit and then slowly approaches the equilibrium 
value, i.e. vanishes. Comparing the different cluster sizes, we can say that the almost 
exact final-state dynamics is obtained with — 10 for times up to imax = 10. This is 
traced back to the fact that due to the spatial separation between U and V vertices, 
vertex corrections decrease with increasing distance and are sufficiently small for 
Lc = 10. Opposed to the CPT, the results for the isolated cluster exhibit a strongly 
oscillatory behavior as a function of time, and finite-size scaling is obviously impossible 
for larger times. 

6. Summary 

For the calculation of dynamical correlation functions by means of Krylov-space 
techniques, it makes a big difi^erence whether the frequency-dependent Green's or 
spectral function is addressed by approximating the resolvent of the Hamiltonian, 
or, on the other hand, the time-dependent correlation function by approximating the 
exponential of the Hamiltonian, followed by a Fourier transformation to frequency 
representation. 

We first summarize the results for the equilibrium (zero-temperature) single-particle 
Green's function. In the first case, the replacement {u — H)~^\io) i->- P{u! — H)~^P\io) 
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represents an approximation that conserves the first n moments of the spectral density 
if P is the projection onto an n-dimensional Krylov space: Namely, the moments are 
related to the coefficients in the high-frequency expansion of the Green's function, and 
this is obtained via the expansion of the resolvent, 1/{uj — H) = YlT=o / ■ Then, 
conservation of the moments results from the fact that, by construction, H^\io) can be 
exactly represented in the Krylov space /C^dio)) if r < n. 

In the second case, the replacement exp{—iHt)\io) i->- P ex-p{—iHt)P\io) is 
numerically exact for a sufficiently short propagation time t. Using several restarts of the 
Krylov construction, this can be exploited to get the time-dependent Green's function 
and finally, by mean of fast Fourier transformation, the frequency-dependent Green's 
function and spectral density. This time-dependent Lanczos algorithm reduces to the 
conventional one if no restart at all is considered. Like the correction-vector method, it 
provides the numerically exact result at the cost of a largely increased numerical effort 
that is roughly proportional to the number of restarts /. A too small /, however, leads 
to spectral densities violating causality. Comparing the conventional with the time- 
dependent Lanczos technique, the most significant deviations are found at the highest 
frequencies in the spectrum since the convergence of ground state and the low-lying 
excited states with increasing n is very fast in the conventional method. Generally, 
the approximation of the resolvent is a severe one if compared with the approximate 
Lanczos determination of the groimd state. On the other hand, due to the substantially 
larger CPU times necessary, one might tolerate this approximation, in particular in the 
context of the DMFT or dynamical cluster-embedding schemes where the error due to 
the finite small number of sites in the effective impurity or cluster model can be more 
severe. 

For the nonequilibrium case, however, an approach based on the approximation of 
resolvents appears as rather ineffective: Since single-particle correlation functions are 
no longer homogeneous in time, Fourier transformation to the frequency representation 
does not help in solving Dyson's equation which is central to any nonequilibrium 
dynamical embedding method. Nevertheless, the contour-ordered Green's function 
can be represented in terms of resolvents, see (21). Treating these by projection onto 
appropriate Krylov spaces which necessarily must be constructed from the excited states 
as initial states, represents a much less controlled approximation compared with the 
equilibrium case. On the other hand, a correction-vector method would provide the 
numerically exact result, but at the cost of the necessity to solve large sparse linear 
systems of equations for each pair of frequencies on a sufficiently dense frequency mesh. 
Apart from that, a three-dimensional Fourier back transformation is required in the 
resolvent-based approach. 

We have therefore suggested a four-step procedure to compute the Green's function 
as a function of t and t' on the Keldysh-Matsubara contour which is numerically exact 
and much more efficient than an approach to resolvents based on correction vectors: 
After (i) finding the initial state |0) of the system as the ground state of an initial 
Hamiltonian by means of a standard Lanczos procedure, (ii) = exp{—iHt)\0) 
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is computed by means of the Krylov technique, followed by (iii) the back propagation 
= exp{iHt)ca\^{t)) for every t and (iv) the evaluation of a scalar product 
{'^a'{t')\'^ait)) from which the different components of the contour-ordered Green's 
function are obtained. Step (iii) is most CPU time consuming. Compared to the time- 
dependent Lanczos approach to the equilibrium Green's function, the main complication 
consists in the fact that the initial state Ca|$(t)) for the time propagation in (iii) is time- 
dependent itself. There are no problems, however, to get all components of Gaa'{t,t') 
for a half- filled Hubbard or single- impurity Anderson model with Lc — 12 sites and 
for t, t' < 100 on the real branches, for example, using a standard PC using the total 
particle number and the ^-component of the total spin as good quantum numbers. 

As a simple application of the nonequilibrium exact-diagonalization solver, we have 
considered the nonequilibrium cluster-perturbation theory [33]. For the single-impurity 
Anderson model in a semi-infinite chain geometry, a magnetic excitation that is localized 
in the vicinity of the correlated impurity is expected to dissipate into the uncorrelated 
and unpolarized bath in the process of time. This is nicely seen within the NE-CPT if 
the hopping V, linking the linear cluster of the first Lc sites with the infinite uncorrelated 
bath, is treated as the inter-cluster hopping by means of all-order perturbation theory 
in V and U. Within the NE-CPT the neglected vertex corrections are controlled by 
the spatial distance between the U and the V vertex and thus by the cluster size. Our 
calculations for Lc = 4 — 10 sites show a systematic improvement and give the essentially 
exact result on a time scale of t < 10 (in units of the inverse hopping). This is just the 
scale which is typically accessible by means of dynamical impurity or cluster embedding 
approaches and which is relevant, for example, to estimate the speed of information 
processing in atomic-scale all-spinbased devices [26] . 

There are several points that may be addressed in future studies: In case of driven 
systems or within the context of nonequilibrium DMFT, where the time dependence of 
the Hamiltonian H{t) is more complicated than a simple sudden quench of a parameter, 
the Krylov time evolution must be carried out using a time discretization step At that is 
considerably shorter than any characteristic time scale of H{t). More efficiently, higher- 
order commutator-free exponential time-propagation algorithms [52] can be applied. 
Another interesting line is the computation of higher-order Green's functions using the 
Krylov approach. Correlation functions depending on four independent time variables 
are required, for example, in the context of the nonequilibrium dual-fermion approach 
[34]. Finally, for the application of the nonequilibrium exact-diagonalization solver 
within the NE-CPT the use of advanced quadrature formulas is promising to extend the 
Keldysh branch, i.e. the limit tmax up to which observables can be traced. Applications 
to two-dimensional lattice models are particularly interesting as there is hardly an 
alternative to an exact-diagonalization solver. 
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